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We present a systematic study of corrections to scaling in the self-organized critical forest-fire 
model. The analysis of the steady-state condition for the density of trees allows us to pinpoint the 
presence of these corrections, which take the form of subdominant exponents modifying the standard 
finite-size scaling form. Applying an extended version of the moment analysis technique, we find 
the scaling region of the model and compute the first non-trivial corrections to scaling. 
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I. INTRODUCTION 

The term self-organized criticality (SOC) [0 refers to 
a set of driven dissipative systems that, under the action 
of a very small external driving, evolve into a critical 
state characterized by avalanches broadly distributed in 
space and time, which lead to divergent (power-law) re- 
sponse functions. Since its introduction by Bak, Tang, 
and Wiesenfeld [||, the concept of SOC has been the 
object of a very intense research activity, covering both 
theory and numerical simulations. 

Among the many models proposed so far exhibiting 
SOC behavior, the forest-fire model (FFM) ||-|] is one 
of the most simply defined and well understood. The 
FFM is a three states cellular automaton defined on a d- 
dimensional hypercubic lattice. Each site of the lattice is 
occupied either by a tree, by a burning tree, or is empty. 
Every time step, the cellular automaton evolves accord- 
ing to the following set of rules: (i) each burning tree 
becomes an empty site; (ii) every tree with at least one 
burning nearest neighbor becomes a burning tree; (hi) a 
tree becomes a burning tree with probability /, irrespec- 
tive of its nearest neighbors; (iv) an empty site becomes a 
tree with probability p. The FFM possesses two charac- 
teristic time scales M : the average time for a tree to grow 
1/p and the average time between fires 1//. In the limit 
of the double infinite time scale separation, 1 > p > /, 
the model displays critical behavior 00]: i.e., fires are 
distributed according to power laws. The magnitudes 
characterizing a fire are the total number of trees burnt 
s, and the total time duration of the fire t (measured as 
the total number of parallel updatings of the algorithm) . 
In the critical state, withp// 3> 1, the probability distri- 
butions of sizes and times have been observed to follow 
the standard finite-size scaling (FSS) hypothesis 0: 

P( S> 9)=s- T 'F{£;), (1) 

W) = ^(^), (2) 

where = p/ / is the critical parameter of the model ||, 
and t x and X x are scaling exponents characterizing the 



critical state Finally, T{z) and Q{z) are cut-off func- 
tions that are constant for z — > and decay exponentially 
fast for z — > 00. 

The precise determination of critical exponents is a rel- 
evant issue in order to firmly establish universality classes 
and the upper critical dimension, that on their turn are 
fundamental in the theoretical understanding of the crit- 
ical nature of the model. While the numerical determi- 
nation of the overall power law behavior is a relatively 
easy task, a very accurate determination of critical expo- 
nents from numerical simulations can suffer from strong 
systematic biases due to the distribution's lower and up- 
per cut-offs. More subtly, the assumption of the FSS 
form does not take into account the presence of correc- 
tions to scaling due to subdominant exponents. These 
corrections are more evident for small values of the var- 
ious magnitudes and for deviations from pure criticality 
(d^ 1 ^ 0). On the other hand, for a sufficiently large 
value of 9, one can safely assume that the scaling (Q)- 
(0) is essentially correct. Let us then define the scaling 
regime of the model by the parameter # S cai, defined such 
that the single scaling picture is correct for > # sca i; 
in principle, 9 SC a.i is a magnitude which depends on the 
microscopic details of the model. However, the value of 
$scai is in general unknown, and when analyzing numer- 
ical data, there is no a priori way to ascertain whether 
the range of 9 at our disposal is large enough. 

In this paper, we will show that in the stationary state 
of the forest fire model, the presence of scaling correc- 
tions arises naturally. The analytical inspection of the 
steady state condition points out the presence of sub- 
dominant scaling corrections and calls for an extended 
scaling framework allowing the evaluation of the scaling 
regime and the various corrections to scaling present in 
the model. The proper treatment of scaling corrections 
permits a more precise estimate of the leading exponents. 
In order to analyze the occurrence of correction to scal- 
ing in a systematic way, we generalize the powerful mo- 
ment analysis introduced in Refs. |H],[llJ to a more gen- 
eral scaling form. Within this framework, we are able to 
estimate the value 9 SC& \, above which the simple form (|l|) 
is meaningful. We thus obtain corrected exponents, and 
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the values of the first subdominant exponents. 

The paper is organized as follows: In Section II, by 
analyzing the steady state condition, we show the in- 
eluctable emergence of subdominant corrections to scal- 
ing in the FFM. In Sec. Ill, we review the moment anal- 
ysis technique, and outline its extension to probability 
distributions with subdominant terms. Sec. IV provides 
numerical evidence of our results by means of extensive 
simulations of the FFM in d = 2. Finally, our conclusions 
are summarized in Sec. V. 



II. STATIONARITY CONDITION AND SCALING 
CORRECTIONS 

The necessity to include corrections to scaling indeed 
arises naturally in the FFM, by just considering the 
steady-state condition of the model . For any value 
of 9, at large times the FFM sets in a steady state char- 
acterized by an average constant density of trees, p t , and 
empty sites, 1 — pt (the density difference after and be- 
fore a fire is negligible, being of order 6/L 2 ). The density 
of trees is known to display the asymptotic behavior at 
large 9 |J| 



Pt = Pt 



a9~ 



(3) 



Computer simulations in d = 2 provide the values pf ~ 
0.408 and a ~ 0.5 ||. In the steady-state, and for 
a fixed 9 value, the average number of growing trees, 
p(l — pt)L d , must equal the average number of burnt 
trees, fp t (s) e L d , where {s) g is the average size of a fire. 
Therefore, the mean number of trees burnt by a fire is 
given by 
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Pt 



(4) 



In the limit 9 — > oo, the mean tree density reaches its 
critical value pf , and we recover the usual relation (s) g ~ 
9 = p/ f ||]. However, for any finite value of 9, the system 
is subcritical and p t < pf . Substituting the expression 
(|J) into (||), we obtain 



<«>6 



9 x 



1 - pf + a9- a 



Pt 



a9- 



PT 



q1 — a 



(pry 



0{9 



l-2a\ 



(5) 



That is, neglecting corrections of order 9 l ~ 2a (which is 
valid since a ~ 0.5), the form of the average fire size for 
finite 9 is 



( s ) e = de + c 2 



(6) 



where the Ck are constants independent of 9. Inspection 
of Eq. (||) proves that it is impossible to obtain such 
a 9 dependence for the average avalanche size with an 
FSS of the form (0). We are therefore forced to admit a 



more complex scaling form. These corrections to scaling, 
which on the other hand are well-known in the field of 
equilibrium and non-equilibrium critical phenomena jfjj , 
take the form of subdominant corrections to the leading 
(infinite 9) scaling form of the probability distributions. 
The most general form of these corrections is 



(7) 



In the last equation, r* and A* are subdominant ex- 
ponents, correcting the infinite 9 behavior, and !F*{z) 
is a cut-off function that decays faster than T(z) when 
z — > oo. In this way, and for fixed 9, the effects of the 
corrections are expected to be more noticeable for small 
values of x. The ellipsis denotes other possible correc- 
tions, which are of lower order compared to the first one. 

In this perspective, a very accurate measurement of 
critical exponents cannot escape the precise knowledge 
of the extent of the intermediate region in which scaling 
corrections are still noticeable. In particular, a method 
of analysis which takes into account the presence of sub- 
dominant exponents is required for a fully consistent 
analysis of the scaling properties at finite values of 9. 



III. MOMENT ANALYSIS 

The determination of the scaling exponents for the 
FFM has been performed most often in previous works 
by a direct measurement of the slope of a log-log plot 
[|],|]]|Jl3|,[l|| • This procedure yields the exponent t s by 
means of a straightforward linear regression. The expo- 
nent A s is then computed by imposing the constraint (|l|) 
for different values of 9, using the previously computed 
value of t s ||. 

Even though with this procedure (sometimes supple- 
mented with extrapolations and/or local slope analysis) 
one can determine the exponents within a 10% accuracy, 
its performance is affected by the existence of the up- 
per and lower cutoffs, which render difficult its applica- 
tion. Moreover, any binning performed to smooth the 
numerical distributions can lead to biases very difficult 
to assess. In this respect, it is better to use analysis tech- 
niques that use the whole set of data (not only the power 
law regime) and contain explicitly the system-size depen- 
dency. In the SOC field, the moment analysis has been 
introduced by De Menech el al. in the context of the 
two dimensional Bak-Tang-Wiesenfeld model [ ]l0"| , [n] ] and 
has been successfully applied to both deterministic and 

In the following we introduce 
extend the method in order to 



stochastic models 15-18 



the moment analysis anc 
deal with scaling forms that make explicit the presence 
of subdominant corrections. 
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A. Single scaling form 

In this section, we concentrate in the moment anal- 
ysis of the fire size, following Refs. |lO| , [lT|| . We start 
with a distribution fulfilling the scaling form (m. The 
g-th moment of the distribution is defined by [s q ) e = 
Jj 00 s q P(s, 9)ds. Inserting the scaling form of P into this 
expression yields the 9 dependence 



( S 1) e = 0A.(*+1-t.) 



>T{z)d 



z, 



(8) 



where we have used the transformation z = s /9 Xb . For 
large values of 9, and provided that q > r s — 1, the lower 
limit of the integral in (||) can be replaced by 0. Wc 
then have (s q ) e ~ Q x s(i+ 1 ~ T s) . In general, we can write 
{s q ) g ~ 9 tJs< - q \ where the exponents <r s (q) can be ob- 
tained as the slope of a log- log plot of {s q ) s as a function 
of 9. Comparing with (||), one has (s q+1 ) g / '(s q ) ~ 9 Xs 
or <r s (q + 1) — <r s (q) = A S) so that the slope of a s (q) as a 
function of q is the cutoff exponent; i.e. A s = da s (q)/dq. 
This is not true for small q, because the integral in m) 
is dominated by its lower cutoff. In particular, the lower 
cutoff becomes important for q < t s — 1. Once the ex- 
ponent A s is known, we can estimate t s from the scaling 
relationship (2 — t s )A s = cr s (l). 

The results of the moment analysis must finally be 
checked by means of a data collapse analysis. The 
initially assumed FSS hypothesis (!]) has to be veri- 
fied, and must be consistent with the calculated expo- 
nents. This can be done by rescaling s — > s/9 Xs and 
P(s,9) -> P(s : 9)9*° T °. Data for different values of 9 
must then collapse onto the same universal curve if the 
FSS hypothesis is to be satisfied. Complete consistency 
between the methods gives the best collapse with the ex- 
ponents obtained by the moment analysis. 



B. Moment analysis with corrections to scaling 

Let us now develop the formalism of the moment anal- 
ysis for a distribution with corrections to scaling of the 
form (^), where we will only keep the first non-trivial 
correction. By plugging this form into the definition of 
the g-th moment, we obtain 



<»*>, 



3 -<+9jr (_£_ ) ds 



= ft \s(q+l-r s ) I z -^+1T{z)dz 



+ 0A;( g +i-r;) / z -<+ir*{z)dz. 



(9) 



In the integrals of the previous expression we have explic- 
itly written the dependence on the lower cut-off. For 9 
sufficiently large and q > max(r s ,T*) — 1, the lower lim- 
its tend to zero, and thus we expect the integrals to be 



independent of 9. However, we cannot discard in general 
a possible dependence on q (through the exponent in the 
integrand). We have therefore 



C{q)9 



\ s (q+l-T s ) 



C*(q)9 



*:(g+i-0 



(10) 



where we have defined the constants (independent of 9) 



C(q) = / z- Ts+q F(z)dz, 
Jo 

POO 

C*(q) = / z- T ' +q T*{z)dz. 
Jo 

Analysis of the general Equation (|l(]) is extremely dif- 
ficult, due to the impossibility to separate the two lead- 
ing behaviors 9 Xsq and 9 X ' q . In order to achieve further 
progress, we must somehow simplify relation ([!(]). To 
do so, we proceed to make an ansatz, whose validity will 
have to be numerically verified a posteriori. The ansatz 
consists in assuming the identity 



A., — A* , 



(11) 



that is, the cut-off exponents do not suffer from correc- 
tions. The physical interpretation of this single cut-off 
exponent for both the leading and subdominant terms 
in the size probability distribution is related to the pres- 
ence of a unique and well-defined divergent characteristic 
size in the avalanche evolution. Under this assumption, 
Eq. (nG) becomes 



C{q)9 



A«(2-r«) 



C*(q)9 



A s (2-<) 



Specializing this relation to q — 1 we obtain 
(s) g = \C(1)9 X ^ 2 - T ^ +C*(l)(9 A ^ 2 -<) 



(12) 



(13) 



Comparing now with the expression for the average fire 
size (first moment), Eq. (|^), we can identify X s (2 — r s ) = 
1 and A s (2 — r*) = 1 — a, from which we obtain the 
exponents 



1/A S 



(1 - a)/X a 



(14) 



Using the previous relations, we can express Eq. ( |l2|) as 
a function of the exponents a and A s alone: 



(A 



- flA s (g-l)+l 



[C(q) + C*(q)9- 



(15) 



Eq. ( |l5| ) suggests the correct strategy to work out the 
moment analysis. Firstly, we observe that the quantity 



r,(0) 



<*«>6 



0A„(g-l)+l 



C{q)+C*{q)9- 



(16) 



depends only on 9~ a . We can use this fact to verify the 
validity of the ansatz A* = A s by plotting T q (9) as a func- 
tion of 9~ a for different values of q, and checking whether 
or not the plots have linear dependence. Secondly, we 
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FIG. 1. Scaling of the average density of trees as a function 
of the parameter 6 — p/f. 
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FIG. 2. Scaling of the average fire size as a function of the 
parameter 6 —p/f. 



note that the second term between brackets in the r.h.s. 
of ([l5|) decreases with increasing 9. For 9 sufficiently 
large, this second term is negligible with respect to the 
constant C(q), {s q ) e has a pure power-law dependence 
and we can proceed to compute X s by means of linear 
regressions. Indeed, this observation allows us to define 
quantitatively the scaling region of the model: assuming 
that the ratio C(q)/C*(q) does not depend strongly on 
q, we define 9 sca \ as the value of the scaling parameter 
for which 



c*(i) , 

£7(1) 



seal 



< r, 



(17) 



with r some (arbitrary) small number. For 9 > sca i, the 
approximation {s q ) g ~ 9 Xs ( q ~ 1 " >+1 C(q), and therefore the 
single scaling form (|l|) is correct, within a precision of 
order r. One can thus proceed to compute the quantity 
Csio) by means of regressions limited to values of 9 > 
6* sca i, determine X s by differentiation and, using Eq. (14), 
estimate the rest of the exponents. 



IV. NUMERICAL RESULTS 

In order to check numerically our arguments, we have 
performed extensive numerical simulations of the FFM in 
d= 2, using the algorithm described in Ref. ||. Starting 
from an arbitrary initial configuration, wc update the lat- 
tice according to the following rules: (i) select at random 
a site in the lattice; if the site contains a tree, burn it 
and all the trees that belong to its same forest cluster; if 
the site is empty, proceed to step ii; (ii) select at random 
9 sites; if a site is empty, grow a tree on it; if it contains 
a tree, do nothing. It is easy to see that these set of rules 
are equivalent to the original definition of the FFM, in 
the limit p = + and finite p/f = 9. For large 9, we 
thus ensure the double infinite time scale separation con- 
dition. The system sizes considered are up to L = 19000 



and the values of 9 range from 128 to 32768. Results are 
averaged over 10 7 nonzero fires. 



A. Average density of trees 

In the first place, we study the average density of trees 
as a function of the parameter 9. After discarding a suf- 
ficiently large number of fires (usually 5 x 10 5 ) to en- 
sure that the system is in a steady-state, we compute the 
average number of trees, per unit area, left after each 
fire. The measured p t is fitted to the functional form 
pf — a9~ a using the Levenberg-Marquardt non-linear 
fitting algorithm [jlSf] . We obtain a critical asymptotic 
density of trees = 0.4084 ± 0.0005, and an exponent 
a = 0.47 ± 0.01, in good agreement with previous results 
Q . In Figure |l| we have checked the asymptotic form of 
the average tree density by plotting \og 10 (pf D — pt) as a 
function of log 10 9. 

In Fig. H we check the validity of Eq. (^J). The parame- 
ters computed, using again a non-linear curve fitting, are 
£7i ~ 0.854, C* 2 ~ 1.973, and a = 0.47 ± 0.01. Again, 
we observe a very good fit to the predicted form. In 
view of this results, we select the value a = 0.47 for the 
computations to follow. 



B. Size probability distribution 

Once we have verified the likelihood of corrections to 
scaling in the first moment of the fire size distribution, 
we proceed to analyze the size probability distribution. 
The first step it to compute the threshold 9 sca \ using the 
criterion (|17|). We arbitrarily fix the parameter r = 0.05; 
for this value, together with the estimates of C\ = C(l) 
and Ci — C*(l) obtained by analyzing (s) g , we estimate 
$scai > 3000. For values of 9 larger than 3000; there- 
fore, the single FSS form (uh can be assumed to be valid, 
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FIG. 3. Plot of cr a (q), computed from linear regressions 
from 9 — 4096 to 32768. The slope yields the exponent 
\ a = 1.09 ±0.01. 

and we can proceed along the standard moment analysis 
technique. In Fig. || we plot the moments <r s (q) com- 
puted from linear regressions of log 10 (s q } e as a function 
of log 10 8, for values of 8 between 4096 and 32768. The 
slope of this plot yields the exponent X s — 1.09 ± 0.01; 
finally using the relations (|lj) with a = 0.47, we obtain 
t s = 1.08 ± 0.01 and r* = 1.51 ± 0.02. A summary of 
results is presented in Table |. 

Once we have computed the exponent A s , we can 
check a posteriori the validity of the ansatz A* = A s . 
We do so by plotting the quantity T q (8)/T q (9 IIL i n ) ~ 
(s q ) g /6» A ^9- 1 )+ 1 as a function of 6>- Q , with a = 0.47, 
for several values of q, Fig. [I| For large values of 9~ a , we 
observe a very good linear relationship. The goodness of 
the fit decreases for large q and large 9 (small 8~ a ) be- 
cause in both cases, the q-th moment is dominated by the 
largest avalanches, of which there is poorer statistics. We 
conclude therefore that indeed the assumption A* = A s 
is indeed well justified for the FFM. 

The presence of corrections of the form (0) make im- 
possible to use a standard data collapse to inspect the 
accuracy of our results in the whole range of 9 values. 
However, for 8 > # SC ai, is reasonable to expect a good 
collapse to the single form (Q). We have plotted this data 
collapse in Figure pL for the integrated size distributions. 
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TABLE I. Critical exponents for the FFM model, obtained 
through the slope analysis, Ref. Q , and by means of the mo- 
ment analysis. Figures in parenthesis denote statistical un- 
certainties. 
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0.01 0.03 0.05 0.07 0.09 
0-0.47 

FIG. 4. Rescaled q-th moment T q (0)/r q (8 m i n ) as function 
of 8~ 0Ar . The good linear fit for small 9 validates the ansatz 
A* = A s . The full lines are guides to the eye. 

The collapse for the exponential tail of the distribution 
is quite remarkable. On the other hand, it is poorer for 
small values of s. This effect is due to the very presence 
of corrections to scaling, whose influence is stronger for 
small s. 

Our method provides values which correct previous es- 
timates (namely, in our notation, t s = 1.14 ±0.03 and 
A s = 1.15 ± 0.03, Ref. §) by a 5%. As explained before, 
the discrepancy is due to the fact that in Ref. || ex- 
ponents were computed by directly measuring the slope 
of the probability distributions, a method which is usu- 
ally less accurate. Our value of A s , on the other hand, 
agrees better with the result reported in Ref. which 
was obtained by a method closer in spirit to the moment 
analysis. 

C. Time probability distribution 

To complete the study of the FFM, we proceed in this 
section to apply the moment analysis to the fire time dis- 
tribution. Here there is no a priori clue about the possi- 
ble existence of corrections to scaling. We will therefore 
assume the simple FSS form (|^) and perform the analysis 
for values of 8 larger than 6* sca i. 

Along the same lines followed for the size distribution, 
we define the q-th time moment (t q ) g = t q P(t,9)dt. 
In this case, we have (t q ) e ~ 9 at ^ q \ with A t = da t (q)/dq 
and r t given by the relation (2— r t )A t = cr t (l). In Figure || 
we plot cr t (q) as a function of q, computed by linear re- 
gression for the largest values of 8. From the slope of this 
plot we obtain At = 0.59 ± 0.01, and using this value on 
the precedent scaling relation, we obtain r t = 1.27±0.01. 
The data collapse with these exponents of the integrated 
time distribution is shown in Figure [H. In this case, and 
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FIG. 5. Data collapse analysis of the integrated fire size 
distribution. 6 = 4096,8192,16384, and 32768. Exponents 
used: A s = 1.09, t s — 1.08. Inset: detail of the tail of the 
distribution. 



FIG. 7. Data collapse analysis of the integrated fire time 
distribution. 6 = 4096,8192,16384, and 32768. Exponents 
used: X t = 0.59, r t = 1.27. 
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FIG. 6. Plot of ot{q), computed from linear regressions 
from 9 — 4096 to 32768. The slope yields the exponent 
Xt =0.59 ±0.01. 

on the contrary to the size distribution, the collapse is 
perfect for all values of t, which proves the irrelevance of 
corrections to scaling in the distribution of this magni- 
tude. 

It is interesting to note that our results match quite 
closely the results in Q, namely At = 0.58 and r t = 
1.27 ± 0.07. This fact is accounted for by the method 
employed in Q to compute At, that is, an analysis of the 
lifetime of the largest fire, T max , as a function of 9. This 
procedure indeed amounts to an estimation of the cut-off 
exponent of the time distribution, and is presumably less 
error-prone that a direct measurement of the initial slope 
of the distribution. 



V. CONCLUSIONS 

In summary, in this paper we have shown that subdom- 
inant scaling corrections are inescapable in the forest-fire 
model. The analytical analysis of the stationarity condi- 
tion shows that scaling corrections to a simple FSS form 
of the fires distribution must be included in order to ac- 
count for the model behavior at finite values of 9. In this 
perspective, we have proposed a method to explore cor- 
rections to the finite-size scaling hypothesis in the forest- 
fire model. The method, based in an extension of the 
moment analysis, allows in principle the determination 
of the scaling regime of the models, as well as the compu- 
tation of the first order corrections to scaling. Applying 
our method, we have been able to compute numerically 
corrected values to the scaling exponents, summarized in 
Table |[ and estimate the first nontrivial corrections. We 
note that our approach is complementary with previous 
studies of deviation from scaling due to finite-size effects 
(small L compared with 9) ]20| . 

As a final remark, it is interesting to point out that the 
present method can also be applied to standard sandpile 
models, defined on a lattice of size L with open boundary 
conditions. In this case, however, the applicability of the 
method is hindered by the availability of a smaller range 
of values of the scaling parameter L. Interestingly, pre- 
liminar results with medium system sizes indicate that 
the ansatz \ s = A* may be violated in sandpiles. This 
fact can be related to the more complex structure of 
avalanches in sandpiles (compared with the percolation- 
like fires in the FFM), that induce the presence of more 
than one characteristic avalanche size. Unfortunately, 
the violation of the ansatz renders the computation of 
the corrections considerably much harder. Work is un- 
derways to explore the full structure of the corrections to 
scaling in sadpiles. 
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